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One year ago I presented a paper at this symposium that described some work we had oeen 
doing with the method of averaging. One of the aspects of that paper had to do with a par- 
ticular choice of nonsingular variables that were called equinoctial variables. In particular, 
we investigated analytical averaging techniques in the equinoctial coordinate frame. Some 
very good results were obtained, and we developed some expressions for fairly low-order 
perturbations. We developed explicit expressions for the disturbing potential for the zonal 
harmonics up to and the third-body harmonics up to (a/Rj , where a is the semimajor 
axis of the orbit and R, is the distance to the third body. 

After that we reconsidered the problem and decided that one of the shortcomings of that 
work was the lack of general expressions for the disturbing potential. This lack was a prob- 
lem for several reasons: First, serious problems were encountered in extending the results 
to higher degree terms. The need to extend the results occurs, for example, in mission anal- 
ysis, where the problem of including the effect of one h^er degree third-body term or in- 
cluding the effect of an additional zonal harmonic occurs frequently. Second, having explicit 
expressions for each term made software development complicated because, for each new 
.'•dditional software was required, and there was the possibility of errors at each sljp. 

These are some of the reasons why it was desirable to develop general results for the disturb- 
ing potential in this coordinate frame, and that is what 1 propose to discuss now. We are 
going to consider the gravitational potential and will use several special functions-the Legen- 
dre and associated Legendre polynomials and the , which are called the derived Legendre 
fun :tions. Recently there has been a fair amount of interest in these functions in the context 
of the gravitational potential. (For example. Pines, 1973.) The various disturbing potential 
expansions are listed below. 


* A significantly expanded version of this paper was presented at the AIAA Aerospace Sciences meeting in January 1975. 
Copies of the preprint are available from the author. 
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Disturbing Function Expansions 
Perturbation nth Term 

Third-body harmonics 


Zonal harmonics 


Arbitrary geopotential 
term (n, m) 


Here we are considering the third-body harmonics and the zonal harmonics, in particular. 

(It was originally planned to have some expressions for the tesseral case, but they are still 
in the process of being checked.) We have introduced the semimajor axis into these expres- 
sions, even though it is not required, because we are going to average these expressions and 
want to use existing results in two-body mechanics. For the third-body harmonics, the 4> 
angle will be the angle between the vector to the satellite and the vector to the third body. 

In the case of zonal harmonics, ^ will be the colatitude. 

In averaging these potentials, a two-step process is used. The purpose of step 1 is to obtain 
expressions for the Legendre functions in terms of equinoctial orbital elements. In step 2, 
the averaging is performed. 

One of the important questions in deriving potentials in terms of classical or equinoctial 
elements is how to get the fast variable motion into the potential. There are several different 
options, and they have a direct effect upon the amount of manipulation and on the com- 
plexity of the final results. Also related to these options is the level of accuracy that can be 
achieved. For third-body harmonics, there are really two options; We can use the mean 
anomaly for both the satellite and for the third body and end up with a very general potential 
This was done by Kaula in 1966. Or we can use direction cosines (relative to an orbital refer- 
ence frame) to get the location of the third body. If direction cosines are used, there is not 
quite the flexibility in modeling that is had with the Kaula approach, but the expressions are 
a lot more compact. In our work with the method of averaging, the use of the direction co- 
sines is quite appropriate. 

For zonal harmonics, we bring in the motion of the fast variable by using the orbital true 
longitude. For the arbitrary geopotential term, by which is meant tesserals, there are iwo 
options that are somewhat related to the two options that exist for third-body harmonics: 

We can use the mean variable for the satellite and Greenwich sidereal time, which is particu- 
larly appropriate for the study of resonance cases; or we can use the true longitude of the 



(C„mCOsmX + S^^sinmX) 
\ = a-d 
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satellite and Greenwich sidereal time, which appears to be most appropriate if we want to 
assume that the central body rotational angle is fixed during the averaging. Such a potential 
might be particularly appropriate for lunar satelUtes. 

Several expressions for cos ^ are listed below. These correspond to various attempts to 
develop the third*body potential. The first expression listed for cos ^ in terms of right 
ascension and declination of the satellite and the third body was used by Kaula. He then 
used the addition theorem and the definition of the inclination function to obtain expres- 
sions for the Legendre polynomials in terms of the classical elements. These results have 
great generality because we can impose all types of resonance constraints on them, yet they 
are also very complex. 

Expressions for Cosine ^ 


cos ill = a cos v + P sin V 

o = ^ P,P = ^ Q 


cos ^ = a cos u 0 sin u 

u = v + cj,o = ^ • N,P = ^ 'M 

cos ^ = a cos L + 0sinL 

L = v + t») + n,a = ^ = - g 

More recent work in third-body perturbations uses the concept of expressing third-body 
motion by direction cosines relative to an orbital coordinate frame. Kaufman in 1970 and 
Lorell and Liu in 1 97 1 , in particular, worked with the apsidal coordinate frame. Subsequently, 
Kozai worked with a coordinate frame where the orbital x-axis was along the nodal aossing 
and the third axis was in the orbit plane. Essentially, he used the argument of latitude as 
his fast variable. At about the same time, we used the true longitude as our fast variable and 
also used the direction cosines of the third body relative to the equinoctial coordinate frame 
(FGW), in which the x-axis points at the origin of the latitudes. 

One of the problems involved in working with Legendre polynomials and, in particular, with 
the direction cosines formulation, was that there did not appear to be an addition theorem 
for that case. In most of the work using that approach, the basic recuraon relationships 
have been used to generate the Legendre functions. There have not been general formulas 
for the required Legendre functions with the argument a cos L + j3 sin L. This is the same 
problem I reported on a year ago. Recently, we took another look at the addition theorem, 
looked at the expressions in terms of the direction cosines, and found that we could reformu- 
late the problem in such a way that the addition theorem still applied; the steps in that deriva- 
tion are shown in equations 1 and 2. 
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First, we introduce a phase angle, L', and deflne the trigonometric functions-cos L* and 
sin L'— to fit the form of the a«iUi>:ion theorem. Then a, 6 , otj, and 63 must be defined. 

a = L cos 5= 1 

03 = L' sin 6=0 ( 1 ) 

00563= \/a^ +0* sin 6 3 = -^ 1 - = 7 

Here those angles are treated as arbitrary quantities. By substituting equation I into the 
addition theorem; 

n 

PJcos^) = PJ0)PJy) + 2 ^ 

m= 1 


(n - m)! 


(n + m)! 


P..(0)Pnn.W“sn’(L-L'). 


( 2 ) 


In equation 2 , 7 is the direction cosine of the third body relative to the vector normal to the 
orbit plane. 

If we define some and polynomials in a and similar to the w • Pines did in 1973, 
the addition theorem can be simplified. The final result is 

n 

P„(cos = P (0) P„( 7 ) + 2 ^ P„J0) q,„(7) ic Jot./o cos m L + SJaJ) sin m Lj (3) 

m= 1 


where 


and 



(4) 

S.c.« = !.(«♦ jor 

(5) 


Because = 0 for n - m equal to an odd integer, gives even powers of the 7 function. 
It can be expressed quite well in terms of a polynomial in a and p. Again, it is finite, so we 
have a general form, a modified addition theorem, that is useful in this case. 

To give an idea of what these polynomials look like, the first one that appears in the third- 
body perturbation is evaluated as follows; 

3 

P-fo cos L + ^ sin L) = ( 37 ^- 1) +— ((a*- cos 2L + 2afi sin 2L] 

2 ' 4 4 

1 r 3 3 1 ( 6 ) 

= y y(o*+/ 3*)+— (o*-U*)cos2L+3o08in2L-ll. 

This checks against the results that were presented in our paper a year ago. It has been 
checked for several cases and gives exactly the results obtained previously through very 
long and arduous manipulations; in fact, some errors were found in those previous results 
through comparison with the general formula. 
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Anotlier major point of this paper is to emphasize the computation of the potential by 
using the recursion formula, and below are listed the recursion formulas (equations 7 to 10) 
tliat are appropriate to the computation of (o cos L + ^ sin L). 


(n+l)P„^,(7) 

= (2n + l)rP^(r)-nP„ ,(7) 

(7) 


(2n - 1) 7 Q„.,,„(7) - (m + n - 1) Q„.j _„(7) 

(8) 


= aCJa,fi).0SJaJ) 

(9) 


= P CJa.fi) 

(10) 


Equation 7 is the standard recursion for the Legendre polynomials, and equation 8 is the 
recursion for the derived Legendre functions (this is given in Ananda and Broucke, 1973, and 
Pines, 1*^73). Equations 9 and 10 are the recursions for the and polynomials, respec- 
tively. 


The final formula for the averaged potential also has the same and polynomials, with 
k and h as the argument, and the same recursions will still apply. We are interested in ob- 
taining averaged potentials, and we want to average with respect to a time-oriented variable, 
so we average with respect to the mean longitude. In equation 1 1 , we take the standard 
formula for the average of (r/a)" cos m v and recast it in terms of the true longitude: 


1 

lit 




e“"‘-dX 


H)”(^ 


(k + jh)"*F 


m -n - 1 


m - n 


2 


,m+ l,h^+k^ 


)• 


( 11 ) 


It should be noted that the formulas are a hypergeometric series and that the polynomial 
previously referred to as appearing in these integrals has the argument k + jh. The symbol 
j is the square root of -1 . 

We wanted to obtain a recursion formula for computing the right-hand sides of these integrals, 
and it happens that one has already been derived and is presented in Cook’s paper on the PROD 
program (1973). The formula for these integrals is as follows: 


i /'( t )'' (- t )' « 


+ ih)'"B|;;^j(h2+k2) (12) 


where 


Tift *»ft*l nic-1 »k (^+k)(8-k-l) + 1^2^ ftk + 1 

=l,Bg -BJ+ (h +k 

4ic (lc+ 1) 
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So along with (a cos L + 13 sin L), the hypogeometric series and the eccentricity poly- 
nomials can be reached recursively. 

The following is a result for the nth term in the third-body potential expressed in terms of 
equinoctial elements and the direction cosines: 



where 


n 

^ * k*) Qn. m (T') fC™ («.0) C™ (k, h) + («,/3) S„ (k, h)l 


m=l 



(n - m)! / n + m + 1 
(n + m)! \ m 


(13) 


All of the terms in equation 1 3 can be obtained recursively. In writing a program, we would 
probably initially compute the coefficients ^ and then store them as data at the beginning 
of the program. We do not obtain all of the terms indicated in the equation. We get approx- 
imately (n/2) terms in that summation because some of the associated Legendre functions 
with zero argument are zero, depending on whether n - m is an even or odd number. We have 
a finite, closed-form expression, and we verified several of the results presented a year ago 
with this formula. The next step is to show how the same analysis is done for zonal har- 
monics. The potential due to the nth zonal harmonic is 

II IRX l>f" 

(t) 

where rl/ is the colatitude. Cos ^ can be expressed in classical elements by 

cos - sin i sin (v a>) . (15) 

In equation 1 5, i is the inclination and v is the true anomaly. We can express cos ^ quite 
well in terms of the equinoctial elements as is shown in the next equation. 

-2p 2q 

cos\l»= cosL+ sin L (16) 

1 + p* + 1 + p* + q^ 

In equation 16, p and q are equinoctial elements. If we assume auxiliary variables defined 
by 
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-2p 


and 




(17) 


2q 

l+p"+q^ 


the argument of the Legendre function in equation 14 has the same form as it did in the 
third-body case, and we can use the addition theorem we derived previously for (o cos L + 
sin L). 


For averaging over the mean longitude, we have to somewhat modify our averaguig formula. 
We could use the previous averaging formula with negative values of n, but that would result 
in an infinite series in the eccentricity, which can be seen by examining the hypogeometric 
series involved. Essentially , we made a quadratic transformation in t he hypogeom etric series, 
changed the arguments a little, introduced a variable x defined as 1/y^l -h^ -k^ , and found 
that we could get a finite series in the eccentricity and x for the required integrals: 


lit 



x*" ‘ (k + ih)™ B™(h^ +k^). 


(18) 


We also use the same polynomials, the polynomials in this case, which can be used in 
computing the potential for both third-body and zonal effects. There will be another saving 
in computation if both effects are included. 


Next we give a general expression for the potential due to a zonal harmonic, that is, the 
averaged potential in terms of the same eccentricity functions, derived Legendre functions, 
and the same a-|3 polynomials; 


a® 


(0)P„(r)x"" ‘ (h^+k^) 


n n 1 ( 

£ *"• O'’ ♦ k’) Q„ „ C) (C„ («,«) C„ (k. h) t <«,» S„ <k. 

m=l 


where 


W. 


n, m 



(n - m)! / n - 1 
(n + m)! \ m 


Again, everything can be derived recursively; terms that are useful in the the third-body case 
are useful here, and terms that are useful for the lower order zonal harmonics can be used to 
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compute the higher order ones. Therefore, I suspect that a very efficient computation pro- 
gram can be built using this formula, with the added advantage of having only one formula 
rather than many to verify. 
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